Moving from Peru to Scotland exposed me to winters that were not only
strikingly beautiful but also harshly challenging. Here is when I first
heard the term ‘Seasonal Depression’, medically classified as Seasonal
Affective Disorder (SAD), being thrown around at this time of the year -
on the cusp of winter. Through casual conversations with friends,
classmates and acquaintances I got to know that almost everyone suspects
they have suffered from SAD at some point throughout their time in
Scotland, whereas they are a local or an international student like
me.
Therefore, this has prompted me to explore whether this widely expressed
sentiment is reflected at a population-level.
SAD is widely discussed in literature, with a focus in populations living at high latitudes, which associates limited winter daylight with increased depressive symptoms (Melrose, 2015). Scotland sits at similar latitudes of well-known Scandinavian cities like Copenhagen, Bergen and Oslo, making this study exciting as studies on Nordic countries have found significant differences in mood changes suggesting a correlation between daylight hours and mood (Adamsson, Laike and Morita, 2018). However, the extent to which these seasonal patterns have been discussed and researched in Scotland is minimal when compared to those on Nordic countries.
This report focuses on investigating whether antidepressant
medication prescribing in Scotland shows seasonal patterns and whether
those patterns are related to regional daylight exposure and deprivation
levels.
The research question asks:
To what extent do seasonal variations in daylight hours affect antidepressant medication prescribing across Scottish NHS Health Boards?
For this investigation, prescription data was taken as a direct inferential measure of the population-level mental health burden.
Null (H0): Antidepressant prescriptions do not have seasonal patterns across NHS Health boards.
Alternative (Ha): Antidepressant prescriptions have seasonal patterns across NHS Healthboards.
Where if Ha is true, the following is proposed:
Prescriptions are higher during seasons with shorter daylight hours and lower during those with longer daylight hours.
This report proposes two complementary hypotheses that further enrich the multi-factorial lens through which seasonal effects influence prescribing:
I.Latitudinal Influence. Health Boards located further north will show stronger seasonal fluctuations in prescribing.
II. Socioeconomic Influence. Health Boards with a higher socioeconomic deprivation (lower SIMD ranking) will have a consistently higher benchmark antidepressant use regardless of season, making the impact of varying daylight hours and temperature be even greater.
Thus, these hypotheses collectively examine how Scottish antidepressant prescription patterns are shaped by seasonal, geographic and social conditions.
Variables were regionally and seasonally categorised according to Met Office UK guidelines. For information on specific categorisation see Appendix 1.1 and Appendix 1.2.
Raw data from multiple institutional websites include all medications prescribed by NHS Health Boards. To differentiate and include only Antidepressant prescriptions, only prescriptions with a BNF item code starting with ‘0403’ were used.
The investigation window used was March 2024 to February 2025
(inclusive of March 1st, 2024 - February 28th, 2025) to assess the
potential cyclical nature of results and make adequate comparisons
between seasons.
All links mentioned below are included in the code within the Methodology.
docs/data. See original TXT files in Appendix 1.3. (Note: Met Office considers Dec
2024, Jan 2025 and Feb 2025 as Winter 2024)docs\data within the repo.library(tidyverse)
library(here) # for the upkeep of the directory structure
library(janitor) # for data cleaning
library(lubridate)
library(gt) # for table building
library(sf) # for geospatial visualisation
library(ggplot2)
library(ggtext)
library(patchwork)
library(plotly)
library(cowplot)
library(stringr)urls_prescr <- list(
prescr_jan_june_2024 <- "https://www.opendata.nhs.scot/dataset/84393984-14e9-4b0d-a797-b288db64d088/resource/f0df380b-3f9b-4536-bb87-569e189b727a/download/hb_pitc2024_01_06-1.csv",
prescr_july_dec_2024 <- "https://www.opendata.nhs.scot/dataset/84393984-14e9-4b0d-a797-b288db64d088/resource/f3b9f2e2-66c0-4310-9b8e-734781d2ed0a/download/hb_pitc2024_07_12-1.csv",
prescr_jan_june_2025 <- "https://www.opendata.nhs.scot/dataset/84393984-14e9-4b0d-a797-b288db64d088/resource/9de908b3-9c28-4cc3-aa32-72350a0579d1/download/hb_pitc2025_01_06.csv")
# Reads all Health Board prescription data in a loop to avoid repetition
prescr_list <- map(urls_prescr,
~read_csv(.x) %>%
clean_names())
# Binds together everything in a single tibble
prescr_raw <- bind_rows(prescr_list, .id = "source_file") %>%
mutate(paid_date_month = str_trim(as.character(paid_date_month))) %>%
select(-source_file)
glimpse(prescr_raw)## Rows: 2,249,380
## Columns: 9
## $ hbt <chr> "S08000015", "S08000015", "S08000015", "S0800001…
## $ dmd_code <dbl> 1.001011e+15, 1.001411e+15, 1.001811e+15, 1.0018…
## $ bnf_item_code <chr> "0603020J0AAAEAE", "1001010P0AAAHAH", "1310012F0…
## $ bnf_item_description <chr> "HYDROCORTISONE 20MG TABLETS", "NAPROXEN 250MG G…
## $ prescribed_type <chr> "VMP", "VMP", "VMP", "VMPP", "VMP", "VMPP", "VMP…
## $ number_of_paid_items <dbl> 25, 53, 275, 1, 181, 2, 487, 1432, 66, 1, 1, 283…
## $ paid_quantity <dbl> 1244, 4046, 4695, 15, 25320, 240, 24924, 65820, …
## $ gross_ingredient_cost <dbl> 145.58, 187.17, 1111.15, 3.55, 4093.40, 38.80, 5…
## $ paid_date_month <chr> "202401", "202401", "202401", "202401", "202401"…Filtering out by “bnf_item_code” to keep only those prescriptions that are antidepressants (codes starting with ‘0403’) and that were prescribed within investigation window.
# Only keeping antidepressant codes and aggregating prescriptions per HB per month
prescr_monthly <- prescr_raw %>%
filter(!is.na(bnf_item_code)) %>%
filter(str_detect(bnf_item_code, "^0403")) %>%
mutate(paid_date_month = as.integer(paid_date_month)) %>%
group_by(hbt, paid_date_month) %>%
summarise(number_of_items = sum(number_of_paid_items, na.rm = TRUE)) %>%
arrange(paid_date_month)
# Subsetting to fit our investigation window
prescr_monthly <- prescr_monthly %>%
filter(paid_date_month >= 202403, paid_date_month <= 202502)
prescr_monthly %>%
summarise(rows = n(), min_month = min(paid_date_month), max_month = max(paid_date_month))## # A tibble: 15 × 4
## hbt rows min_month max_month
## <chr> <int> <int> <int>
## 1 S08000015 12 202403 202502
## 2 S08000016 12 202403 202502
## 3 S08000017 12 202403 202502
## 4 S08000019 12 202403 202502
## 5 S08000020 12 202403 202502
## 6 S08000022 12 202403 202502
## 7 S08000024 12 202403 202502
## 8 S08000025 12 202403 202502
## 9 S08000026 12 202403 202502
## 10 S08000028 12 202403 202502
## 11 S08000029 12 202403 202502
## 12 S08000030 12 202403 202502
## 13 S08000031 12 202403 202502
## 14 S08000032 12 202403 202502
## 15 SB0806 1 202412 202412spr_months_202425 <- c(202403, 202404, 202405)
sum_months_202425 <- c(202406, 202407, 202408)
aut_months_202425 <- c(202409, 202410, 202411)
win_months_202425 <- c(202412, 202501, 202502)
seasons_202425 <- tibble(
paid_date_month = c(spr_months_202425, sum_months_202425, aut_months_202425, win_months_202425),
season = c(rep("Spring", length(spr_months_202425)),
rep("Summer", length(sum_months_202425)),
rep("Autumn", length(aut_months_202425)),
rep("Winter", length(win_months_202425))))# Health Board official NHS names list
hb_names <- read_csv("https://www.opendata.nhs.scot/dataset/9f942fdb-e59e-44f5-b534-d6e17229cc7b/resource/652ff726-e676-4a20-abda-435b98dd7bdc/download/hb14_hb19.csv") %>%
clean_names()
# Met Office regional mapping
north_hb <- c("NHS Highland", "NHS Western Isles", "NHS Orkney", "NHS Shetland")
east_hb <- c("NHS Borders", "NHS Lothian", "NHS Fife", "NHS Tayside", "NHS Grampian", "NHS Forth Valley")
west_hb <- c("NHS Ayrshire and Arran", "NHS Dumfries and Galloway", "NHS Greater Glasgow and Clyde", "NHS Lanarkshire")
metoffice_hb_region <- tibble(
hb_name = c(north_hb, east_hb, west_hb),
region = c(rep("North", length(north_hb)),
rep("East", length(east_hb)),
rep("West", length(west_hb))))
# Health Boards per region according to Met Office scottish territorial classifications
hb_regional <- hb_names %>%
full_join(metoffice_hb_region, by = "hb_name") %>%
select(-c(hb_date_archived, hb_date_archived, hb_date_enacted, country))From October 2024 data file. The data file chosen was deliberate as October 2024 marks approximately half-way of the investigation window.
hb_pop <- read_csv("https://www.opendata.nhs.scot/dataset/e3300e98-cdd2-4f4e-a24e-06ee14fcc66c/resource/cec9341e-ccba-4c71-afe4-a614f5e97b9f/download/practice_listsizes_oct2024-open-data.csv") %>%
clean_names() %>%
select(hb, sex, all_ages) %>%
filter(!sex %in% c("Male", "Female")) %>%
group_by(hb) %>%
summarise(hb_population = sum(all_ages, na.rm = TRUE)) %>%
ungroup()prescr_seasonal <- prescr_monthly %>%
full_join(hb_regional %>%
select(hb, hb_name, region), by = join_by(hbt == hb)) %>%
full_join(hb_pop, by = join_by(hbt == hb)) %>%
full_join(seasons_202425, by = join_by(paid_date_month))full_join() there were 4 rows with NA in
prescription data and population. Upon further investigation and a look
at the object hb_names, these were shown to have had their
hbt numbers archived in 2018 and 2019. Thus, they were removed from the
dataset.prescr_seasonal <- prescr_seasonal %>%
filter(!is.na(hb_name)) %>%
filter(!is.na(paid_date_month))
# checking if there is a missing region or population
prescr_seasonal %>%
filter(is.na(region) | is.na(hb_population)) # tibble of 0 x 7 shows we have eliminated all NAs## # A tibble: 0 × 7
## # Groups: hbt [0]
## # ℹ 7 variables: hbt <chr>, paid_date_month <dbl>, number_of_items <dbl>,
## # hb_name <chr>, region <chr>, hb_population <dbl>, season <chr>A calculation of items_per_1000_people was made to allow for population weighting of prescription items. This is because all NHS Health Boards have different populations, thus, comparing their “number_of_items” solely would have affected by population sizes.
prescr_seasonal_standard <- prescr_seasonal %>%
mutate(items_per_1000 = (number_of_items/hb_population)*1000)This is somewhat challenging. Met Office sunshine datasets are
provided as text files (.txt) with white spaces as delimiters. I
converted them into .csv files using Excel to ensure reproducibility and
easier data handling. These files can be found in the
docs\data folder in the repo. A custom function was written
to: a) extract the 2024 values for each region, b) reshape seasonal
columns into long format, and c) standardise season naming. This
produced a tidy dataset of seasonal total daylight
hours per region. > Note: these .csv files contain columns named
spr, sum, aut, win
and a year column for each year and season
# Reading all CSV files first
daylight_north <- read_table(here("docs", "data", "R_north_scotland_sunshine.csv")) %>%
clean_names()
daylight_east <- read_table(here("docs","data", "R_east_scotland_sunshine.csv")) %>%
clean_names()
daylight_west <- read_table(here("docs","data","R_west_scotland_sunshine.csv")) %>%
clean_names()
# Making a function to avoid code repetition for each .csv file
daylight_season_function <- function(data, region_name, year_filter, season_cols = c("win", "spr", "sum", "aut"), year_cols = c("year_12", "year_13", "year_14", "year_15"), full_season_names = c("Winter", "Spring", "Summer", "Autumn")) {
# built-in checker
if(length(season_cols)!= length(year_cols)) stop("season_cols and year_cols need to have the same length")
if(length(season_cols)!= length(full_season_names)) stop("full_season_names and season_cols need to have the same length")
# processing each individual season
season_list <- map2(season_cols, year_cols, ~ {
data %>%
select(all_of(.x), all_of(.y)) %>%
filter(.data[[.y]] == year_filter) %>%
rename(year = all_of(.y))})
# joining all four seasons together
season_complete <- reduce(season_list, full_join, by = "year") %>%
relocate(all_of(season_cols), .after = last_col()) %>%
mutate(across(all_of(season_cols), as.numeric)) %>%
pivot_longer(cols = all_of(season_cols), names_to = "season", values_to = "daylight_hrs") %>%
mutate(season = recode(season, !!!set_names(full_season_names, season_cols)),
region = region_name) %>%
arrange(year, factor(season, levels = full_season_names)) %>%
filter(!is.na(daylight_hrs))
return(season_complete)}daylight_season_functiondaylight_season_function() was used to wrangle and
select specific data from each Met Office regional daylight dataset.
all_seasons_daylight contains all total daylight hours per
season per region for the duration of the investigation window.
# Using the function for each region and making sure the "year" category is a character
season_daylight_north <- daylight_season_function(daylight_north, "North", "2024")
season_daylight_north <- season_daylight_north %>%
mutate(year = as.character(year))
season_daylight_east <- daylight_season_function(daylight_east, "East", "2024")
season_daylight_east <- season_daylight_east %>%
mutate(year = as.character(year))
season_daylight_west <- daylight_season_function(daylight_west, "West", "2024")
season_daylight_west <- season_daylight_west %>%
mutate(year = as.character(year))
# Merging all regional daylight data for each season
all_seasons_daylight <- bind_rows(season_daylight_north, season_daylight_east, season_daylight_west)
all_seasons_daylight## # A tibble: 12 × 4
## year season daylight_hrs region
## <chr> <chr> <dbl> <chr>
## 1 2024 Winter 122. North
## 2 2024 Spring 375. North
## 3 2024 Summer 330. North
## 4 2024 Autumn 230. North
## 5 2024 Winter 164. East
## 6 2024 Spring 362. East
## 7 2024 Summer 453. East
## 8 2024 Autumn 269. East
## 9 2024 Winter 137. West
## 10 2024 Spring 367. West
## 11 2024 Summer 403. West
## 12 2024 Autumn 230. West# Median rank per Health Board was used to summarise the distribution of deprivation
simd_raw <- read_csv("https://www.opendata.nhs.scot/dataset/78d41fa9-1a62-4f7b-9edb-3e8522a93378/resource/acade396-8430-4b34-895a-b3e757fa346e/download/simd2020v2_22062020.csv") %>%
clean_names()
simd_hb <- simd_raw %>%
select(hb, simd2020v2rank) %>%
group_by(hb) %>%
summarise(SIMD_median_rank = median(simd2020v2rank, na.rm = TRUE)) %>%
ungroup()
# Joining daylight, SIMD and seasonal prescription data
analysis_prescr_simd <- prescr_seasonal_standard %>%
full_join(all_seasons_daylight, by = join_by(season, region)) %>%
full_join(simd_hb, by = join_by(hbt == hb))
# Checking if any NAs are present
analysis_prescr_simd %>%
summarise(missing_daylight = sum(is.na(daylight_hrs)), missing_simd = sum(is.na(SIMD_median_rank))) # result should show a 14 x 3 tibble with value zero (0) for each## # A tibble: 14 × 3
## hbt missing_daylight missing_simd
## <chr> <int> <int>
## 1 S08000015 0 0
## 2 S08000016 0 0
## 3 S08000017 0 0
## 4 S08000019 0 0
## 5 S08000020 0 0
## 6 S08000022 0 0
## 7 S08000024 0 0
## 8 S08000025 0 0
## 9 S08000026 0 0
## 10 S08000028 0 0
## 11 S08000029 0 0
## 12 S08000030 0 0
## 13 S08000031 0 0
## 14 S08000032 0 0hb_shp_geo <- st_read(here("docs","data", "NHS_healthboards_2019.shp")) %>%
clean_names()## Reading layer `NHS_healthboards_2019' from data source
## `/Users/florenciasolorzano/Documents/data_science/B218332/docs/data/NHS_healthboards_2019.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 14 features and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 7564.996 ymin: 530635.8 xmax: 468754.8 ymax: 1218625
## Projected CRS: OSGB36 / British National Gridanalysis_geo_prescr <- analysis_prescr_simd %>%
group_by(hbt, season, SIMD_median_rank) %>%
summarise(av_items_per_1000 = mean(items_per_1000, na.rm = TRUE), av_daylight = mean(daylight_hrs)) %>%
ungroup() %>%
full_join(hb_shp_geo, by = join_by(hbt == hb_code)) %>%
st_as_sf()Figure 1: A seasonal summary table by region
full_seasonal_table <- analysis_prescr_simd %>%
group_by(region, season) %>%
summarise(av_daylight = mean(daylight_hrs, na.rm = TRUE), av_items_per_1000 = mean(items_per_1000, na.rm = TRUE)) %>%
ungroup() %>%
arrange(region, factor(season, levels = c("Spring", "Summer", "Autumn", "Winter")))
full_seasonal_table %>%
mutate(av_items_per_1000 = round(av_items_per_1000, 2),
av_daylight = round(av_daylight, 2)) %>%
gt(groupname_col = "region") %>%
cols_label(
season = md("Season"),
av_daylight = md("Mean Total Daylight (hrs)"),
av_items_per_1000 = md("Mean Prescriptions (units/1000 people)")) %>%
tab_header(
title = md("Antidepressant Prescriptions per 1000 and Total Daylight hours by region"),
subtitle = "March 1st, 2024 - February 28th, 2025") %>%
fmt_number(columns = c(av_items_per_1000, av_daylight), decimals = 2)| Antidepressant Prescriptions per 1000 and Total Daylight hours by region | ||
| March 1st, 2024 - February 28th, 2025 | ||
| Season | Mean Total Daylight (hrs) | Mean Prescriptions (units/1000 people) |
|---|---|---|
| East | ||
| Spring | 361.80 | 113.06 |
| Summer | 452.60 | 113.84 |
| Autumn | 268.60 | 113.79 |
| Winter | 163.60 | 114.51 |
| North | ||
| Spring | 374.60 | 121.67 |
| Summer | 330.40 | 120.00 |
| Autumn | 229.70 | 121.45 |
| Winter | 121.90 | 121.96 |
| West | ||
| Spring | 367.40 | 137.25 |
| Summer | 403.10 | 138.44 |
| Autumn | 230.30 | 138.47 |
| Winter | 137.20 | 139.04 |
Figure 1 shows the mean seasonal total daylight hours and the mean antidepressant prescriptions per 1000 population for each region.
Figure 2 shows a dual bar chart graph displaying antidepressant prescriptions per 1000 people proportional to each NHS Health Board population and average daylight hours per season, faceted by Scottish Geographical Region (North, East, and West).
dual_bar_analysis <- analysis_prescr_simd %>%
group_by(region, season) %>%
summarise(av_items_per_1000 = mean(items_per_1000, na.rm = TRUE), av_daylight = mean(daylight_hrs, na.rm = TRUE)) %>%
ungroup() %>%
mutate(season = factor(season, levels = c("Spring", "Summer", "Autumn", "Winter"))) %>%
pivot_longer(cols = c(av_daylight, av_items_per_1000),
names_to = "variable",
values_to = "value")
dual_bar_plot <- dual_bar_analysis %>%
ggplot(aes(x = season, y = value, fill = variable, text = ifelse(
variable == "av_daylight",
paste0("Daylight (hrs): ", round(value, 2)),
paste0("Prescriptions: ", round(value, 2))))) +
geom_col(position = position_dodge(width = 0.9), alpha = 0.8) +
facet_wrap(~region, nrow = 1, scales = "free_x") +
scale_y_continuous(
name = str_wrap("Average Total Daylight (hrs)", width = 30),
breaks = seq(0, max(dual_bar_analysis$value), 50),
sec.axis = sec_axis(~ ., name = str_wrap("Average Antidepressant Prescriptions (units/1000 people)", width = 35), breaks = seq(0, max(dual_bar_analysis$value), 50))) +
scale_fill_manual(values = c("orange", "skyblue"),
labels = c("Daylight (hrs)", "Prescriptions (units/1000)")) +
labs(
title = str_wrap("Average seasonal total daylight hours and average antidepressant prescription items by region", width = 63),
subtitle = "Prescriptions per 1000 people across Scottish Health Boards",
x = "Season",
fill = "") +
theme_minimal(base_size = 13) +
theme(axis.text.x = element_text(angle = 45, hjust = 1),
legend.position = "top",
legend.title = element_blank(),
legend.text = element_text(size = 10),
strip.background = element_rect(fill = "gray90", color = NA),
strip.text = element_text(face = "bold", size = 12),
plot.title = element_text(face = "bold", size = 16, margin = margin(b = 5), hjust = 0.3),
plot.title.position = "panel",
plot.subtitle = element_text(face = "bold", size = 12, margin = margin(b = 10), hjust = 0.5),
panel.grid.major.x = element_blank(),
panel.grid.minor = element_blank(),
plot.margin = margin(t = 20, r = 10, b = 40, l = 10, unit = "pt"))
dual_bar_plot
Prescriptions seem to only marginally increase from Spring to Winter across all regions. Interestingly, summer across regions don’t show the lowest average antidepressant prescriptions nor does winter the highest prescriptions. Rather, it seems that antidepressant prescription trends increase as the seasonal year progresses.
Figure 3 shows geospatial visualisation of prescriptions per season per Health Board, where two choropleth maps have been faceted to facilitate seasonal comparisons.
map_prescr_seasons <- analysis_geo_prescr %>%
ggplot() +
geom_sf(aes(fill = av_items_per_1000), size = 0.15, color = "darkgrey") +
scale_fill_distiller(palette = "Blues", direction = 1) +
facet_wrap(~season, nrow = 1) +
labs(title = "Seasonal Antidepressant Prescriptions (March 2024 - February 2025)", subtitle = "Prescriptions by Scottish Health Board per 1000 people comparable with Average Total Daylight (hrs)", fill = "units/1000 people") +
theme_void() +
theme(
plot.title = element_text(size = 14, face = "bold", hjust = 0.5),
plot.subtitle = element_text(size = 10, margin = margin(t = 10, b = 20), hjust = 0.5),
legend.title = element_text(face = "bold", size = 6))
map_daylight_seasons <- analysis_geo_prescr %>%
ggplot() +
geom_sf(aes(fill = av_daylight), size = 0.15, colour = "darkgrey") +
scale_fill_distiller(palette = "Oranges", direction = 1) +
facet_wrap(~season, nrow = 1) +
labs(fill = "Av. Total Daylight (hrs)") +
theme_void() +
theme(legend.title = element_text(face = "bold", size = 6))
full_map_plot <- map_prescr_seasons / map_daylight_seasons +
plot_layout(heights = c(1,1))
full_map_plot
As seen, there are season daylight changes specific to each region
with Western Scotland having the lowest throughout all seasons. It seems
like the southernmost Health Boards are situated, the higher the
anitdepressant prescriptions per 1000 people there are. The choropleth
map shows this almost uniformly despite seasonal changes with varying
daylight hours. The northernmost regions are not the ones with most
prescriptions despite having the lowest total average daylight hours
overall. However, Shetland and the Western Isles show to have higher
prescription units than the Highlands.
Figure 4 shows geospatial visualisation of whether SIMD rankings and antidepressant prescriptions are correlated.
# Bivariate map bins
biv_bins <- analysis_geo_prescr %>%
mutate(prescr_bin = ntile(av_items_per_1000,3),
simd_bin = ntile(SIMD_median_rank, 3),
biv_class = paste0(prescr_bin, "-", simd_bin))
# bivariate map bin colours
biv_palette <- c(
"1-1" = "#e8e8e8",
"2-1" = "#b8d6be",
"3-1" = "#64acbe",
"1-2" = "#d4b9da",
"2-2" = "#a5add3",
"3-2" = "#4a7bb7",
"1-3" = "#c994c7",
"2-3" = "#df65b0",
"3-3" = "#dd1c77")
# bivariate map bin matrix
biv_matrix <- expand.grid(prescr_bin = 1:3, simd_bin = 1:3) %>%
mutate(simd_bin = simd_bin, biv_class = paste0(prescr_bin, "-", simd_bin))
# bivariate map legend
biv_legend <- biv_matrix %>%
ggplot(aes(x = prescr_bin, y = simd_bin, fill = biv_class)) +
geom_tile(color = "white") +
scale_fill_manual(values = biv_palette, guide = "none") +
scale_y_continuous(breaks = 1:3, labels = c("Low", "", "High")) +
scale_x_continuous(breaks = 1:3, labels = c("Low", "", "High")) +
labs(x = "Prescriptions", y = "SIMD Rank") +
coord_fixed(ratio = 1)+
theme_minimal(base_size = 9) +
theme_void()+
theme(
axis.title = element_text(size = 8, face = "bold"),
axis.text = element_text(size = 6),
panel.grid = element_blank(),
plot.margin = margin(t = 0, r = 5, b = 0, l = 0))
# bivariate map
map_biv_solo <- biv_bins %>%
ggplot() +
geom_sf(aes(fill = biv_class), size = 0.1, colour = "white") +
scale_fill_manual(values = biv_palette, guide = "none") +
facet_wrap(~season, nrow = 1) +
labs(
title = str_wrap("Seasonal Antidepressant prescriptions comparison with Median Deprivation Rankings", width = 63),
subtitle = str_wrap("Prescriptions per 1000 people across Scottish NHS Health Boards | Scottish Index of Multiple Deprivation (SIMD) Rank"),
fill = "Bivariate classification") +
coord_sf(expand = FALSE) +
theme_void()+
theme(
plot.title = element_text(face = "bold", size = 14, margin = margin(r = 0.5, b = 0.5), hjust = 0.5),
plot.title.position = "panel",
plot.subtitle = element_text(size = 10, margin = margin(t = 10, b = 20), hjust = 1),
panel.spacing = unit(0.1, "lines"),
strip.text = element_text(face = "bold", size = 8),
plot.margin = margin(t = 0, r = 10, b = 0, l = 10, unit = "pt"))
full_bivariate_map <- plot_grid(map_biv_solo, biv_legend,
ncol = 2,
rel_widths = c(4, 1),
align = "t")
full_bivariate_map
Lower SIMD ranks indicate more deprived Health Board populations. The general trend seems to be that Health Boards that have a higher SIMD rank have consistently low prescribing (NHS Highland, Grampian). Eastern Scotland shows the most deprivation and high prescriptions. Hence, there is an indication of that higher baseline prescribing is delivered more in deprived areas during the season with least sunshine throughout the year. This is also corroborated in the Figure 5 in Appendix 1.4. A thing to note is that NHS Shetland, despite having a high SMID Rank shows high prescriptions, which remains unchanged seasonally. Which perhaps suggests there might be a seasonal influence that was not shown through the analysis done previously.
There is marginal evidence suggesting that antidepressant prescriptions increase during the winter period in several Health Board regions across Scotland in Figure 1. This remains an inference as statistical testing must be done to confirm significant differences between seasons. Figure 3, nonetheless, showcases health board differences better but without any added prescription differences. All of which prompt the acceptance of the null hypothesis given that there is not enough evidence showing concrete differences between seasonal prescriptions in Figure 2 and Figure 3.
Furthermore, Health board regional (latitudinal) differences exist to a limited degree. Northernmost boards with the least total average daylight show do not show striking seasonal patterns, even when comparing summer and winter data. These are not visually consistent with complementary hypothesis I, which must be not accepted.
Finally, deprived areas (lower SIMD rank) somewhat show higher baseline prescribing and marginal seasonal changes in Figure 4. However, the introduction of this factor demonstrate that Deprivation superimposes seasonality, suggesting the social factors are also equally if not more important than environmental factors when researching depressive disorders within the population as a whole. And hence, that the complementary hypothesis II cannot be accepted without further research..
Limited studies in Scotland suggest that 9% of the population suffers SAD yearly but they themselves recognise that more research and quantitative visibility must be done (Boddon et al., 2024). Hence, the conclusions of this report must not be interpreted as the non-existence of SAD in Scotland but rather, as a reflection of the condition’s complexity in visibility and awareness.
Adamsson, M., Laike, T. and Morita, T. (2018). Seasonal Variation in Bright Daylight Exposure, Mood and Behavior among a Group of Office Workers in Sweden. Journal of Circadian Rhythms, 16(1).. doi:https://doi.org/10.5334/jcr.153.Melrose, S. (2015).
Melrose, S. (2015). Seasonal Affective Disorder: An Overview of Assessment and Treatment Approaches. Depression Research and Treatment, [online] 2015(1), pp.1–6. doi:https://doi.org/10.1155/2015/178564.
Boddon, S., Lorimer, H., Parr, H. and Williams, C. (2024). SAD geographies: Making light matter. Progress in geography, 48(5). doi:https://doi.org/10.1177/03091325241252846.
# Appendix 1.1. Season Categorisation as per Met Office UK Guidelines
# Spring: March, April, May
# Summer: June, July, August
# Autumn: September, October, November
# Winter: December, January, February
# Appendix 1.2. NHS Health Board Geographical Categorisation
#This categorisation was made according to the Met Office's territorial delineation of Scotland based on the distribution of climate measurements as available in their website.
# Northern Scotland: NHS Highland, NHS Western Isles, NHS Orkney, NHS Shetland
# Eastern Scotland: NHS Borders, NHS Lothian, NHS Fife, NHS Tayside, NHS Grampian, NHS Forth Valley
# Western Scotland: NHS Ayrshire and Arran, NHS Dumfries and Galloway, NHS Greater Glasgow and Clyde, NHS Lanarkshire
# Appendix 1.3. Met Office Seasonal Data (TXT files)
# Northern Scotland:
# <https://www.metoffice.gov.uk/pub/data/weather/uk/climate/datasets/Sunshine/ranked/Scotland_N.txt>
# Eastern Scotland:
# <https://www.metoffice.gov.uk/pub/data/weather/uk/climate/datasets/Sunshine/ranked/Scotland_E.txt>
# Western Scotland:
# <https://www.metoffice.gov.uk/pub/data/weather/uk/climate/datasets/Sunshine/ranked/Scotland_W.txt>Figure 5 - Interactive scatter plot
scatter_plot_data <- analysis_geo_prescr %>%
st_drop_geometry() %>%
filter(!is.na(SIMD_median_rank), !is.na(av_items_per_1000))
scatter_plot <- plot_ly(
data = scatter_plot_data,
x = ~SIMD_median_rank,
y = ~av_items_per_1000,
type = "scatter",
mode = "markers",
color = ~hb_name,
text = ~paste0(
"HB:NHS ", hb_name, "<br>",
"Season: ", season, "<br>",
"Prescriptions/1000: ", round(av_items_per_1000,1), "<br>",
"Av. Total Daylight (hrs): ", round(av_daylight, 1)),
hoverinfo = "text",
marker = list(size = 10, opacity = 0.8)) %>%
layout(
title = list(text = "Antidepressant Prescriptions according to Median Deprivation Ranks",
font = list(size = 14),
x = 0.5),
annotations = list(list(
x = 0.5,
y = 1.05,
text = "Prescriptions per 1000 people across Scottish NHS Health Boards | Scottish Index of Multiple Deprivation (SIMD) Rank",
xref = "paper",
yref = "paper",
showarrow = FALSE,
font = list(size = 10),
xanchor = "center",
yanchor = "bottom")),
xaxis = list(title = "SIMD Rank (Higher = Less Deprived)"),
yaxis = list(title = "Antidepressant Prescriptions (units/1000)"),
margin = list (t = 90, r = 100, b = 85, l = 85))
scatter_plot